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ABSTRACT 

The perturbations of weakly-viscous, barotropic, non-self-gravitating, New- 
tonian rotating fluids are analyzed via a single partial differential equation. The 
results are then used to find an expression for the viscosity-induced normal-mode 
complex eigenfrequency shift, with respect to the case of adiabatic perturbations. 
However, the effects of viscosity are assumed to have been incorporated in the un- 
perturbed (equilibrium) model. This paper is an extension of the normal-mode 
formalism developed by Ipser & Lindblom for adiabatic pulsations of purely- 
rotating perfect fluids. The formulas derived are readily applicable to the per- 
turbations of thin and thick accretion disks. We provide explicit expressions 
for thin disks, employing results from previous relativistic analyses of adiabatic 
normal modes of oscillation. In this case, we find that viscosity causes the funda- 
mental p- and g-modes to grow while the fundamental c-mode could have either 
sign of the damping rate. 

Subject headings: hydrodynamics — accretion, accretion disks — stars: oscilla- 
tions 
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1. INTRODUCTION 

The problem of perturbations of rotating fluids is of considerable interest in astrophysics. 
These perturbations can give us information not only about stars and accretion disks (e.g., 
Friedman & Ipser 1992; Nowak & Wagoner 1992; Kato et al. 1998), but also about the 
properties of the compact object at the center of the latter (e.g., Perez et al. 1997). 

Traditionally, this problem has been approached in two different ways. The classical 
way has been to use a Lagrangian displacement formalism (Lynden-Bell & Ostriker 1967). 
An alternative method uses two scalar potentials (proportional to the perturbations of the 
pressure and gravitational potential) within the framework of Eulerian fluctuations (Poincare 
1885; Ipser & Managan 1985; Managan 1985). In the spirit of the latter approach, Ipser 
& Lindblom (1991b, hereafter IL; 1992) have developed an elegant formalism to analyze 
adiabatic pulsations of perfect Newtonian (as well as relativistic) fluids that are stationary, 
axisymmetric, and purely rotating (i.e., have no meridional component of velocity). They 
look for normal-mode solutions and reduce the problem to two coupled second-order partial 
differential equations for the two potentials. 

The inclusion of viscosity in such a system is an important step in the modeling of 
various types of astrophysical phenomena (e. g., Ipser & Lindblom 1991a). The purpose of 
this paper is to explore some of the consequences arising from the presence of weak viscosity, 
taking into consideration that its effect is two-fold: a) it changes the value of the unperturbed 
fluid variables that characterize the equilibrium configuration, and b) it changes the form of 
the differential equations that govern the evolution of the perturbations. However, we assume 
that step (a) has already been carried out. We will work within the Cowling approximation 
(i.e., neglecting gravitational perturbations), which reduces the system to a single differential 
equation. This approximation is often appropriate (Cowling 1941; Ipser & Lindblom 1992), 
and is certainly valid for most accretion disks. 

In section 2 we present a derivation of the governing partial differential equation. In 
section 3 we use this result to obtain a formula for the viscosity-induced shift of the complex 
eigenfrequencies of the perturbations. This enables us, in section 4, to compare our results 
with those of Ipser & Lindblom (1991a), who perform the calculation for the simpler case 
of pure uniform rotation of the equilibrium model (which then can contain no effects of 
viscosity). In section 5, we apply our formalism to thin accretion disks. Explicit integral 
expressions (involving the eigenfunctions) are obtained for the contributions to the growth 
(or damping) rate of a mode. We then evaluate the order of magnitude of that rate for the 
fundamental modes. Except where indicated, we shall employ the notation and conventions 
of IL. 
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2. GENERALIZED IPSER LINDBLOM FORMALISM 

The equations of motion of a Newtonian fluid are 

d tP + V a ((w a ) = 0, (2-1) 
p(d t v a + v h V b v a ) = - V a p + pV a $ + V b a ab . (2-2) 

In the cases we shall consider, $ represents a specified fixed, axisymmetric gravitational 
potential (which does not necessarily obey the Newtonian field equation). The viscous stress 
tensor a ab is defined by 

In this initial investigation, we adopt a generalized barotropic equation of state which speci- 
fies the pressure p as well as the shear and bulk viscosity coefficients 7] and ( as functions of 
the mass density p only. In these equations, d t and V a represent the partial derivative with 
respect to time and the spatial Euclidean covariant derivative. Below, d a = A = d/dx a in- 
dicate spatial partial derivatives, with x a = r, 0, z. The Euclidean metric g ab and its inverse 
g ab are used to raise and lower tensor indices. 

We assume that the viscosity coefficients are small enough that the effects of the viscous 
term in equation (2-2) can be treated using standard perturbation techniques. To be explicit, 
we assume that the term Vf,a ab has a relative magnitude of order a C 1 compared to the 
dominant terms in the equation. This implies that the viscosity coefficients have a magnitude 
of order ap*LlQ*, where p*, L*, and f2* refer to the typical density, length, and frequency 
scales of the system. (For an accretion disk, corresponds to the thickness of the disk.) In 
this way we can treat a as a constant perturbation parameter with respect to the inviscid 
(a = 0) case, and work to first order in a. 

The presence of viscosity will cause a change in the equilibrium values of the dynamical 
variables. In particular, it can induce meridional components (r, z) of the velocity field. Using 
the noncoordinate basis of orthonormal cylindrical coordinates (with unit vectors r a , <p a , and 
z a ; and = 5^), the latter can be expressed as 

v a = rVL(j) a + v r r a + v z z a . (2-4) 

The last two terms on the right-hand-side of this equation are assumed to have a relative 
magnitude of order a with respect to the first term. (For instance, this is true for a thin 
accretion disk.) 

The equations for the evolution of the Eulerian perturbations (in the Cowling approxi- 
mation) are 

d t Sp + v a V a 5p + 5pV a v a + V a { P 5v a ) = (2-5) 
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and 



<9 t <5t; a + w b V fe (5w a + 5v b V b v a = L + JL—L - JL.1 — 

p p z p z 

+ 



+ 



Specifying the functions 



p- 1 {v b [r ] (V a 5v b + V b 5v a )] + V b [5r](V a v b + VV)] 
V a [(C - lv)V c 5v c ] + V a [(5( - l5v)V c v c }} . (2-6) 

r(r,z)^* (2-7) 
p ap 

f^,*) = -j- (2-8) 
?7 ap 

effectively closes the system of differential equations. (It turns out that the contributions of 
5( are of higher order in a.) The inviscid case has been studied in an elegant fashion by IL in 
more generality. (They included gravitational perturbations and allowed a general equation 
of state, but neglected meridional velocities.) The purpose of this paper is to explore the 
effects of viscosity on their results under the restrictions of the Cowling approximation and 
our generalized barotropic equations of state. 

We look for normal mode solutions to the equations, which (because of the stationarity 
and axisymmetry of the equilibrium configuration) have a time dependence of the form 
e tcrt and an azimuthal angular dependence of the form e tm< ^. The constant a is the inertial 
frequency of the mode and the axial mode number m is an integer. Let us express equation 
(2-6) in a way that shows which part depends explicitly on a and which does not: 

[L\ + a AL a b )5v b = {R a + a AR a )5p . (2-9) 

The derivation of the a = terms, R a and L a b , is found in IL. Here we merely quote their 
results. The operator R a is equal to — V a p _1 , where throughout this paper the operator V a 
is to affect all the terms appearing to its right. The operator L a b is related to the operator 
Q a b of IL by L a b = i(Q~ l ) a b , where (in our orthonormal cylindrical coordinates r, 0, z) 

' u 2 -2iVLu r(n 2 ) yZ 

iuK 2 /2VL uj 2 irtutt iZ (2-10) 
to 2 — k 2 



Q 



Uj{uJ 2 — K 2 ) 



The corotating frequency is uj = a + mf2(r, z), and k 2 = 2f2(2f2 + rQ tr ) is the square of the 
radial epicyclic frequency. We note that uj{uj 2 — k 2 ) is the determinant of Q~ x . (Also note 
that IL use a cylindrical coordinate basis and that their roles of a and uj are switched with 
respect to ours.) The derivation of the expressions for aAL a b and aAR a is straightforward 
and can be found in the Appendix. 
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The key to the procedure is the fact that Q is purely algebraic. This allows us to solve 
for the velocity perturbation: 

Sv c = (L~ 1 ) c a (R a + aAR a )5p - a (L~ 1 ) c a AL a b 5v b . (2-11) 
Thus, to the same order in a, 

5v c » {L~ 1 f a {R a + a Ai2°)«$p ~ a (L-^AL^L' 1 )' d R d 5p . (2-12) 
Terms of order a 2 and higher have been dropped (and will be ignored from now on). 

Taking into account equation (2-7), we can express equation (2-5) in the following way: 

(u - iaC)(p/ P T)5p - iV a (p5v a ) = , (2-13) 

where the operator 

aC = r-\rv r ) yr + u* + v r d r + v z d z . (2-14) 

Throughout this paper, the operators d r (= d/dr) and d z (= d/dz) are to affect all the 
terms appearing to their right. Substituting equation (2-12) in equation (2-13), we obtain 

V a (pQ a 6 V h p-H P ) + ^Sp = a [v a (pQ a b AR%) - t\/ a ( P Q a b AL b c Q c d V d P^p) + tC^dp . 

P P (2-15) 

We have thus reduced the original problem to a single partial differential equation for one 
unknown, 5p. Once solved, 5p, 5r), and 5v a can be obtained via equations (2-7), (2-8), and 
(2-12), respectively. Our master equation (2-15) plays the same role as equation (22) of IL, 
which has the same structure within the Cowling approximation. 



3. EIGENFREQUENCY SHIFTS 

Following IL, for reasons that will become apparent below we will use a new variable: 

5V = 5p/puj . (3-1) 
With this definition we can express equation (2-15) as 

(C + aC v )5V = , (3-2) 

where 

C = W a pQ\W b uj + u 2 p 2 / P T , (3-3) 

C v ee -W a pQ a b AR b up + tV aP Q a b AL b c Q c d V d uj - tCup 2 / P T . (3-4) 
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The operators £ and £ v depend upon a implicitly through the (r and z dependent) properties 
v a , p,p, T), ( of the unperturbed model as well as through the eigenfrequency a. However, we 
assume that the equilibrium model already incorporates the effects of viscosity, at least 
through first order in a. Thus we shall be concerned only with the effects of viscosity on the 
perturbations, which then changes the operator £ only through its dependence on a. 

Expanding the eigenfrequency a = a + Act + . . . and the eigenfunction 5V = 5V + 
aSVi + . . ., we can express equation (3-2) in the following way, to first order in a: 

[£ + (d£/da)Aa + a£ v ] ao (5V + aSV^ = . (3-5) 

All operators are now evaluated at the adiabatic eigenfrequency a . Using the fact that 
C5V = 0, we then obtain 

a C5V 1 + Aa(dC/da)5V + a C V 5V = . (3-6) 

It turns out that with the boundary condition p = on a smooth, compact, finite area 
surface of the body, C is Hermitean (IL). This implies that 

J 5V * £ 6Vi d 3 x = J (£ 5V )* 5V l( i 3 x = 0. (3-7) 

Thus if we multiply equation (3-6) by 5Vq and integrate, we obtain 

Aa(d£/da) + a (£ v ) = ; (3-8) 

where 

(0)= J 5V*O5V d 3 x , (3-9) 
only for operators O evaluated at a = a . 

Then we find (from the structure of £ and £ v ) that 

Re(Aa) = , ilm(Aa) = ~^^y ■ (3-10) 

Equation (3-10) allows us to compute the viscosity-induced shift of the eigenfrequencies, 
provided we are given an equilibrium model. Equation (3-6) can then be used to obtain the 
eigenfunction correction 8V\. 



4. PERTURBATIONS OF UNIFORMLY ROTATING BODIES 



In order to compare our results with those obtained by Ipser & Lindblom (1991a), we 
need to invoke their assumptions of uniform rotation and vanishing of meridional velocity 
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(i.e., we need to set f2, r = Q tZ = v r = v z = 0). Thus, there are no effects of viscosity on 
their equilibrium model since the viscous shear tensor vanishes. In addition, since we do not 
consider self-gravitational effects, we need to set 5$ = in their formalism. Under these 
assumptions, equation (3-10) reproduces their results [equations (27), (28), and (29)] for 
the imaginary part of the frequency induced by viscous effects. [Their r _1 corresponds to 
our Im(Ao-), and their canonical energy E, to our quantity (cu/2)(dC/da).] When dealing 
with differentially rotating bodies, our approach is more straightforward than an extension 
of theirs, which employs a Lagrangian approach. 



5. APPLICATION TO ACCRETION DISKS 



The results obtained so far can be applied to the special case of a thin accretion disk. Its 
physical properties allow us to retain relatively few terms from the ones contained in equation 
(3-10). Adiabatic oscillations of such accretion disks have been studied within a radial WKB 
approximation (Nowak & Wagoner 1992; Perez et al. 1997; Silbergleit & Wagoner 1999) 
A r <r t , with A r > for most of the low-lying modes. The length scales 2h* and r* are the 
typical thickness and radial size of the disk region of interest, and A r is the radial scale of 
variation of the eigenfunction. In addition, we shall only consider values of the axial mode 
number \m\ < 1. 

The system is characterized by the small parameter e = ft,*/r*. In order to use our 
equations we will also need a model for the unperturbed disk, which we take to be similar 
to the one developed by Kita & Kluzniak (1998). In agreement with their results, we take 
v r ~ ae 2 rQ and v z ~ ev r . 

Taking into account these considerations, equation (3-10) gives (to lowest order in a 
and e) the viscous damping rate 



fF„d- 



x 



(5-1) 



To obtain this expression, we employed (many) integrations by parts, with the boundary 
conditions that 5V , rj, and ( vanish. The integrands are 



P 



2up \ J^\SV \ 2 + k 2 S 2 
Tp 



dr 



F 1 = U0[u0 2 {27 1 + i) + K 2 7]} 

F 2 = oo-\2r 1 + 



d f s d5V 
dr V dr 



d 2 5V n 



dz 2 



(5-2) 

(5-3) 
(5-4) 
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^3 

F 4 



u[n(u;- 2 + S 2 (uj 2 + k 2 )) + 2S( V + 0] 



d 2 5V 



drdz 



—ujS 



rr/p ( <9f2 2 \ d 2 ^ 1 
Tp \ dr J dz 2 



d5V 



dr 



(5-5) 
(5-6) 



We have introduced S = l/{oo 2 — k 2 ) and £ = ( 
(3-3), which reduces to the form 



\rj. The IL operator defined by equation 



d 



c Uj2 ^ >2 2 

Tp ^'dr 



S 



d_ 

dr 



+ 



d 



P- 



d 



dz V dz 



(5-7) 



has also been employed. It is interesting to note that the terms proportional to v r and v z do 
not appear in this limit, since their contribution is of higher order in e. This result justifies 
their neglect in previous analyses of thin disks. 

Let us now consider in more detail the effect of viscosity on observationally relevant 
fundamental modes of relativistic accretion disks (near a black hole or neutron star). Since we 
are employing a Newtonian analysis, but use the properties of the relativistic eigenfunctions 
5Vo, our results are only approximate. The trapping of the g- and p-modes is produced by 
the fact that relativistic effects make the radial epicyclic frequency n less than the angular 
velocity of free-particle circular orbits. This reduction produces an inner edge of the disk 
where re(r) vanishes, with the maximum value of re(r) at a slightly larger radius. The trapping 
of the c-modes is related to the reduction of the vertical epicyclic frequency below the angular 
velocity, produced by the angular momentum of the central mass. 



• The g-modes are trapped where the comoving frequency \u\ < k, near the maximum 
of K,(r). They are mathematically similar to the internal gravity (buoyancy) modes of 
stars. For the fundamental (m = 0) g-mode (Perez et al. 1997), d r ~ e 1 l 2 d z , d z ~ 
and 8$ = when applied to 5V . Single terms in F and F 3 dominate, giving 



r = 2 J P K 2 S 2 \d r 5V \ 2 d 3 x 



with S ~ l/(eo" 2 ). The sign of this result agrees with that obtained by Nowak & 
Wagoner (1992), although the magnitude of their growth rate was greater. 

The (inner) p-modes are trapped between the inner edge of the disk and the increasing 
portion of «(r), with > k. They are similar to the pressure (acoustic) modes of 
stars. For the fundamental (m = 0) p-mode (Nowak & Wagoner 1991), d r ~ e^ 3 d z , 
d z ~ 1/h*, and d<f, = when applied to 5Vq. The numerator is dominated by the same 
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term in F 3 as above, as well as F 2 , giving 

1~ J\{2 V + Ov~ 2 \d 2 z SVo\ 2 + v{v 2 + K 2 )S 2 \d r d z 5Vo\ 2 }d 3 x aSl, 



t 2 f[pc- 2 \SV \ 2 + pn 2 S 2 \d r SV \ 2 }d 3 x (ff/fi*) 2 + e 2 / 3 ' 

(5-9) 

with c 2 = Tp/p and S ~ 1/a 2 . The first (second) integral in both the numerator and 
the denominator dominates when e is smaller (greater) than (a/Q*) 3 <C 1. The sign of 
this result agrees with that obtained in a local analysis by Kato et al. (1998), although 
its form differs somewhat from theirs. 

• The c-modes are trapped between the inner edge of the disk and the radius where 
the Lense-Thirring frequency is the eigenfrequency. For the fundamental (m = ±1) 
c-mode (Silbergleit & Wagoner 1999), d<f, — ±1, d z ~ but d 2 z = when applied to 
5Vq. Thus this mode is vertically incompressible. One finds that 

1 f(F 3 + F 4 )d 3 x {K\ 2 „ 

- : '» — <>. . (5-10) 



r 2 J pc- 2 Lu\5V \ 2 d 3 x \X 
where S ~ oj~ 2 . Note that this mode may either grow or be damped by viscosity. 

We remind the reader that the real part of the corrections to all eigenfrequencies vanish. 
Future generalizations of this work should include the effects of buoyancy (nonbarotropic 
equation of state) and effects of higher order in e = /i„/r t . 

This research was supported in part by NASA grants NAG 5-3102 to R.V.W. and NAS 
8-39225 to Gravity Probe B. We thank John Friedman, Lee Lindblom, and Alex Silbergleit 
for helpful suggestions. 



A. APPENDIX: Determination of AR a and AL a b 

The components of V b a ab in orthonormal cylindrical coordinates are 

V b a rb = 2(^, r + r-M7 ? (r-y^ + ^-r-V)], + [r ? (^ + t;f r )], z 

+ 2 W ~\v r >r - r"V - r-V ) + {(C - lv)[r~\rv r ), r + r~ l v^ + v* z ]}, r ; (Ala) 

V b a* b = [ri(v* - r~V + r-yjlr + 2r)r~\v^ r - r~V + r" V ) + 2 r - 2 [ V (v r + 

+ [ V (v* + r-'v^U + r-H(C - h)[r-\rv^), r + r~\\ + ; (Alb) 

V b * zb = fo(t£ + v* r )],r + W-\v\ z + <) + 2(wl), + vr-\v« + r-y^U 

+ {(C - h) [r-\rv r ), r + r~V, + „«]},, . (Ale) 
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These expressions and their perturbations allow us to compute the operator C v . There are 
two terms in the right- hand- side of equation (2-6) that contribute to a AR. We will consider 
each in order. The first contribution comes from —p~ 2 5pV b cr ab . The order a part of this 
contribution can be obtained from the above equations, noticing that — when acting 
on an equilibrium quantity and that the variables v r , v z , r], (, 5r), and 5( are all of order a: 



aAR A = — - 





(rp/ P r)(r]Q !Z ) iZ + (d r r + 2)( V pQ r /pr) 




(A2) 



The second contribution comes from p- 1 V b [5rj(V a v b + VV)]: 



aAR B 



imripQ r /pT 
rd z ( W n z /pT) + (d r r + 2)( W n r /pT) 
imr]pQ tZ /pT 



(A3) 



The term p 1 V a [(5C ~ l^)V c f c ] is of higher order in a and can therefore be ignored. Thus 



a AR = a AR A + a AR B . 



(A4) 



For the determination of a AL, let us start with the contribution from the left-hand-side 
of equation (2-6). The components of the latter are 



d t Sv r + v b V b 5v r + 5v b V b v r = ia5v r + imtt5v r - 2fi<fo* + v r Sv r r + v z Sv r z 



+ v r M + v r z Sv z , 



(A5a) 



d t 5v^ + v b V b 5v' 1, + 5v b V b v' t) 
d t 5v z + v b V b 5v z + 5v b V b v z 



iaSv* + imttdv* + Q5v r + v*5v r + vHv z 

,i ,z 

+ v r 5vt + v'Sv* + r" V<?u* , (A5b) 



ia5v 2 + imVt5v z + v r 5v z + v z 5v z + v z 5v r 
+ v z 5v z . 



The first contribution to a AL comes from the terms with v r or v z : 

Yr V\ z 

a ALa = (v r d r + v z d z ) I + v r /r 

v z v z 



(A5c) 



(A6) 



where I is the identity matrix. Let us turn now to the contributions to a AL coming from 
the right-hand-side of equation (2-6). The first of such contributions comes from the term 
p- 1 V b [r](V a 5v b + V b 5v a )\: 

'2d r r]d r + 2r\d r r~ x — m 2 r~ 2 rj + d z i]d z imr\r~ x {d r — 3r _1 ) 

a AL B = - (d r + Ar^^mr/r^ 1 (d r r + 2)r]d r r~ 1 — 2m 2 r~ 2 r] + d z r/d z 

(d r + r^rjdz imi]r~ 1 d z 
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d z i]d r 
imr~ l d z ri 
(d r + r~ l )rjd r + 2d z r]d z — m 2 r~ 2 rj 



(A7) 




(A8) 



a AL = a ALa + ol AL b + a AL C ■ 



(A9) 
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